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Abstract 

Geometric integration of non-autonomous classical engineering prob- 
lems, such as rotor dynamics, is investigated. It is shown, both nu- 
merically and by backward error analysis, that geometric (structure pre- 
serving) integration algorithms are superior to conventional Runge-Kutta 
methods. 

Key-words: Geometric numerical integration, splitting methods, ro- 
tor dynamics 

1 Introduction 

In this paper we study geometric numerical integration algorithms for non- 
autonomous systems. By classifying the appropriate Lie sub-algebra of vector 
fields, the standard framework for backward error analysis can be used to explain 
the superior qualitative behaviour of geometric methods based on the splitting 
approach. 

The current section continues with a brief review of the general framework for 
geometric methods, mainly following the approach by Reich [2]. In Section[2]we 
study geometric integration of linear systems with non-constant periodic coeffi- 
cients. A numerical example from classical rotor dynamics is given. Conclusions 
are given in Section [3j 

We adopt the following notation. & denotes a phase space manifold of 
dimension n, with local coordinates x = (xi, . . . , x n ). In the case when 
is a linear space we also use P. Further, X(^) denotes the linear space of 
vector fields on 3P. The flow of X € is denoted ip^, where t is the 

time parameter. The Lie derivative along X is denoted Lx- If X, Y € X(3 g ) 
then the vector field commutator [X, Y]x = &xY supplies X(&) with an infinite 
dimensional Lie algebra structure. Its corresponding Lie group is the set Diff (3 g ) 
of diffcomorphisms on with composition as group operation. 

Remark 1.1. More precisely, the group Diff(5 a ) has the structure of a Frechet 
Lie group. See [TJ [3] for issues concerning infinite dimensional Lie groups. 
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As usual, the general linear group of n x n-matrices is denoted GL(n) and 
its corresponding Lie algebra fll(n). We use [A, B] 3 i for the matrix commuta- 
tor AB - BA. 

If V is a metric linear space, then the linear space of smooth periodic func- 
tions R — > V with period 2tt/H is denoted Cq(V). Notice that this space is closed 
under differentiation, i.e., if / G Csi(V) then it also holds that /' G Csi(V). 

1.1 Geometric Integration and Backward Error Analysis 

Let Xs{^) be a sub-algebra of X(3 g ), i.e., a linear sub-space which is closed 
under the commutator. Its corresponding sub-group of Diff(3 g ) is denoted 
Diffs(^). Let X G Xs(^) be a vector field which is to be integrated nu- 
merically. Assume that X can be splitted as a sum of explicitly integrable 
vector field also belonging to X s {&). That is, X = Y + Z where Y, Z G X s {&>) 
and (pyi^z can be computed explicitly. By various compositions, various nu- 
merical integration schemes for Lp l x are obtained. The most classical example 
is $/i = (p'y 2 ° f'z ° <y5y^ 2 , which yields a second order symmetric method (h is 
the step-size parameter of the method). Since (fiyj'Pz £ Diffs(<^), and since 
Difi C 5(J 2> ) is closed under composition (since it is the group operation), it holds 
that $/j G Diffsi^ 1 )- Thus, the splitting approach yields structure preserving 
methods, which is a key property. 

Backward error analysis for structure preserving integrators deals with the 
question of finding a modified vector field X G %s{&) such that $^ = ip%. In 
conjunction with perturbation theory, such an analysis can be used to study 
the dynamical properties of For splitting methods, backward error analysis 
is particularly simple as the modified vector field, at least formally, is obtained 
from the Baker-Campbell-Hausdorff (BCH) formula. For details on this frame- 
work we refer to Reich [5]. 

2 Linear Systems 

In this section we study non- autonomous systems on a linear phase space P with 
global coordinates x = (a?i, . . . , x n ). More precisely, let G be a Lie sub-group 
of GL(n) and g its corresponding Lie sub-algebra. We consider systems of the 
form 

x = A(t)x + f(t) (1) 

where A G Cn(fl) and / G Co(P) is a smooth vector valued periodic function 
with period T — 2n/Vl. Our objective is to construct geometric integrators 
for (jT). Of course, since the system is linear, there is a closed form formula for 
its solution. However, in engineering applications, e.g. finite clement analysis, 
the system is typically very large so computing the exponential matrix, which is 
necessary for the exact solution, is not computationally efficient. Also, it might 
not be possible to analytically integrate / and A over t, which is necessary for 
the exact solution. 

In order to study dynamical systems of the form (JlJ in the framework of 
geometric integration, we need, first of all, to extend the phase space to P = PxM 
to include the time variable in the dynamics. Coordinates on P are now given by 
(x, t) and the new independent variable is denoted r (in practice we always have 
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t(r) = t). Further, we need to find a Lie sub-algebra of X(P) which captures 
the form (|T|). For this purpose, consider the set of vector field on P given by 

£o(P, 0) = {X £ X(P) | X(x, t) = (A(t)x+f(t), a), A e e n (fl), / e C (P), a e 

(2) 

We now continue with some results concerning properties of £,q(P,q). The 
first result states that it actually is a Lie sub-algebra. 

Proposition 2.1. The set of vector fields £q(P,q) is a Lie sub-algebra o/X(P). 

Proof. We need to check that £n(P,g) is closed under vector operations and 
under the Lie bracket. That is, X, Y £ £n(P,{j) should imply aX + bY £ 
£n(P,£|) for a,b £ R and [X,Y] X £ £n(P,fl). 

With = (A(t)cc + f(t),a) and F(x,£) = (B(t)x + g(t),/3) we get 

(aX + bY)(x,t) = ((aA + bB)(t)x + (af + bg)(t),aa + b/3) which is of the desired 
form. Further, 



^ Yh = U(t) A'(t)x + f(t)\ (B{t)x + g(t) 



B{t) B'(t)x + g'(t)\ (Ax + f(t) 
Ma 



(A(t)B(t) - B(t)A(t) + pA'(t) - aB'{t))x + A(t)g{t) - B(t)f(t) + pf(t) - ag'{t) 



which is of the desired form since AB-BA+/3A' -aB' = [A, B^+pA' +aB' £ 
e n (g)a,nd(Ag-Bf + pf-ag')ee Q {P). □ 

From the proof above we immediately obtain the following corollary. 

Corollary 2.1. The set to = Cn(fj) x Cji(P) x R equipped with the induced 
vector operation 

a(A,f,a) + b(B,g,(3) = (a A + bB, af + bg, aa + 6/3), a, 6 £ R 

and with the bracket operation 

[(A, /, a), (B,g, /3)] £ = ([A, B] GL + (3 A' + aB',Ag - Bf + /3/' - 0) 

«s a Lie algebra which is isomorphic to £q(P, q) with isomorphism Iq 9 (A, f, a) i— > 
(AB + /,a)€£n(P,fl). 

Since Cq(P) and Csi(fl) are infinite dimensional it follows that Iq, and there- 
fore also £o(P,g), is infinite dimensional. However, a finite dimensional sub- 
space of Cn(P) is given by 

k 

e n ,*(P) = {/ € e n (P) | /(*) = a + J2 a * cos(iflt) + b, sm(iQt), a,, 6, e P} 

i=i 

(3) 

which is the sub-space of Cq(P) with angular frequencies bounded by kil. Notice 
that the dimension of Cn,fc(P) is (2k + l)n and that Cq !OC (P) = Gq(P) and 
Co,o(P) = P- Further, Cn,fc(P) is closed under differentiation. Clearly, these 
results also holds for the corresponding sub-space Cn,j(fl) of Cn(fl), except that 
the dimension is given by (21 + 1) dim g instead. 

By replacing Gq(P) with (3n,fc(P) and Cn(fj) with Gq,i(q) we get the sub- 
spaces [ n ,fc,/ = Cn,i(fl) x e n , m (P) x R of [ n . In general fo,fc,i is no£ a sub-algebra, 
due to the fact that A,B £ Gq.i(q) does not in general imply AB £ Cn.z(fl)- 
However, in some special cases the implication holds true. 
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Proposition 2.2. The sub-spaces l n ^ kfi = gx Gn,fc(P) xR and [o,fe : oo = Cn(fl) x 
Cn,fe(P) x K are Lie sub-algebras o/[q. Further, ln,k,a is finite dimensional with 
dimension dimg + (2fc + l)n + 1. 

Clearly, ln,fc,o an d (n,fc,oo induces corresponding Lie sub-algebras £n,fc,o(Pjfl) 
and £n,fe,oo(P, 0) of £ n (P, g). 



2.1 Geometric Integration 

In this section we describe an approach for geometric integration of of systems 
of the form ([lj. The approach is based on splitting. To this extent we write 
as an extended system 



^_ /x\ = fA(t)x + f(t) 



X(x,t). (4) 



It is clear that X G £n(P,fl)- Since £o(P,g) is a Lie sub-algebra of X(P) it 
corresponds to a Lie sub-group Diffo of Diff(P). Geometric integration of Q 
now means to find a one-step integration algorithm $^ £ Diffn. There are 
of several ways to obtain such integrators. One of the simplest, but yet most 
powerful ways, is to use a splitting approach. That is, to split the vector field X 
as a sum of two vector fields each of them of the form El . 
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Value 


Unit 


m 


1 


kg 


k 


1 


N/m 


n 


1.02 


rad/s 


e 


0.1 


m • kg 


Xq 


(0,0,0,0) 


(m,m,m/s,m/s) 



Table 1: Data used in the simulations of the rotor dynamical problem. 



2.2 Example: Linear Rotor Dynamics 



This example is the simplest possible rotor dynamical 
problem. It consists of a disc attached to a shaft which 
is rotating with constant angular velocity ft. The shaft 
is held by a bearing, which is modelled as a linear spring 
with stiffness k. (See figure.) The disc is slightly unbal- 
anced, i.e., its centre of mass does not align with rota- 
tional axis. This implies a time-dependent periodic cen- 
trifugal force acting on the rotor. 

The phase space for this system is given by P = K 4 , 
with coordinates x — {qi,q2,Pi,P2), which is the horizon- 
tal and vertical position of the shaft in a plane perpen- 
dicular to the axis of rotation, and their corresponding 
momenta. The equations of motion are of the form Q with 

/ m' 1 \ 
i n ™ - 1 




A = 







and 



/(*) = e& 






- cos(fii) 
sin(ftt) J 



where m is the total mass and e is the magnitude of the unbalance. 

It holds that A T J + J A = 0, so A is an element in the canonical symplectic 
Lie sub-algebra of g((4), i.e., we have q — Sp(4). Further, since A is independent 
of t, and / only contains a single frequency, the appropriate Lie sub-algebra 
of X(E 4 ) is £o i i ! o(R 4 ,sp(4)), which is finite dimensional. 

The eigenvalues of A are ±iy/k/m. Thus if f2 is close to a multiple of the 
eigen frequency to = \fk~jrn of the system starts to resonate. In this example 
we investigate how well various numerical integrators capture that behaviour, 
both qualitatively and quantitatively. 

For the data given in Table |2.2| the problem is numerically integrated with 
four different methods, two which are geometric and two which are not. 



Method 



Geometric? 



Implicit midpoint Yes 

Splitting method Yes 

Heun's method No 

Implicit extrapolation method No 

The results of the xi-variable are shown in Figure [T] Notice that the geo- 
metric integrators captures the resonance phenomena in a qualitatively correct 
way, whereas the non-geometric methods does not show the correct behaviour. 
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Figure 1: Results for qi variable for the simple linear rotor example. From 
top: (1) exact solution. (2) implicit midpoint, geometric. (3) Stormer-Vcrlet, 
geometric, (4) implicit Runge-Kutta, non-geometric, (5) explicit Runge-Kutta, 
non-geometric. All methods are second order accurate. Notice the superior 
qualitative and quantitative behaviour of the geometric methods. 
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3 Conclusions 



A structural analysis of non-autonomous systems has been carried out using 
the framework of Lie sub-algebras of the Lie algebra of vector fields. As a 
direct application, backward error analysis results are obtained for this class 
of problems. Numerical examples of a classical rotor dynamical problem show 
that the geometric methods preserving the structure of the problem indeed are 
favourable over non-geometric dito. 
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